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We develop a microscopic theory for reaction-difusion (R-D) processes based on a generahzation 
of Einstein's master equation with a reactive term and we show how the mean field formulation 
leads to a generalized R-D equation with non-classical solutions. For the n-th order annihilation 
reaction A + A + A + ... + A ^ 0, we obtain a nonlinear reaction-diffusion equation for which we 
discuss scaling and non-scaling formulations. We find steady states with either solutions exhibiting 
long range power law behavior showing the relative dominance of sub-diffusion over reaction effects 
in constrained systems, or conversely solutions with finite support of the concentration distribution 
describing situations where diffusion is slow and extinction is fast. Theoretical results are compared 
with experimental data for morphogen gradient formation. 
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I. INTRODUCTION 



The random walk is the classical paradigm for the microscopic mechanism underlying diffusive processes as demon- 
strated in 1905 by Einstein who showed how the diffusion equation follows from the mean field formulation of the 
microscopic random walk. Here we generalize the formulation for situations where the diffusing particles are also 
subjected to a reactive process. From the phenomenological viewpoint, when diffusion and reaction are coupled, 
these processes are described by reaction-diffusion (R-D) equations. For instance, the evanescence process {A — 0) 
of suspended particles diffusing in a non-reactive medium the concentration of species A, c{r;t), is described by the 
classical R-D equation 

^^c{r;t)=D-^c{r;t) - kc{r;t) , (1) 

where D denotes the diffusion coefficient and k the reaction (evanescence) rate. This classical equation yields a steady 
state solution showing spatial exponential decay of the concentration (but one can equally consider the distribution 



function): c (r) = c(0) exp y — \J k/ D |r| j when particles are injected with a constant flux at r = 0. 

However there are many systems observed in nature where it seems logical to use the language of reaction-diffusion, 
but where non-classical distributions are found, i.e. the steady state spatial distributions are non-exponential e.g. 
when the particles encounter obstacles or are retarded in their diffusive motion, or because the reactive process is 
hindered or enhanced by concentration effects. Such situations are ubiquitous in chemical, rheological, biological, ... 
systems - a typical example being the diffusion and degradation of a morphogen in cells during the early developing 
stage [l| - and are certainly as commonly observed as those that can be described by the idealized R-D system of 
Eq.([T]). This is why approaches to a more general description of R-D phenomena have been proposed, and recent 
developments in this direction [2, 0| are based (i) on a generalization of the diffusive mechanism accounting for time 
delay effects or obstacles hindrance using the continuous time random walk (CTRW) model and corresponding to 
a fractional Fokker-Planck equation (FFP) or the fractional Brownian motion (FBM), and (ii) on a space and time 
dependence of the reaction rate {k — >■ k{r;t)). However the resulting expressions for the steady state distribution 
have so far been subject to controversial comments expressing that "CTRW theory is compatible with available 
experiment" (3] and "that fractional Brownian motion is the underlying process" or that "experimental results 
cannot be explained by a continuous time random walk" Q and "exclude fractional Brownian motion as a valid 
description" [7[. So the present state of the art certainly appears somewhat confusing while it seems nevertheless 
clear that a general R-D theory requires a generalization for both diffusion and reaction. 

Here we present an alternative approach by developing a microscopic theory generalizing Einstein's master equation 
with a reactive term and we show how the mean field formulation leads to the nonlinear R-D equation with non- 
classical solutions. For the n-th order annihilation reaction A + A + A + ... + A 0, we obtain the nonlinear 
reaction-diffusion equation (with no drift) 

^^c{r;t) = ^D-^c"ir;t)-kc-{r;t) , (2) 

for which we discuss scaling and non-scaling formulations and the corresponding range of values of the nonlinear 
exponents. We obtain steady state solutions of the form c(r) = c(0) (1 + Ca.n{D, k) r/u) ^ where v = i^zr^, giving 
long range power law behavior (for n > a) showing the relative dominance of sub-diffusion over reaction effects in 
constrained systems, or conversely (for n < a < 7i + 1) leading to finite support of the concentration c(r) describing 
the situation where diffusion is slow and extinction is fast. An experimental example of morphogen gradient formation 
is discussed. 



II. GENERALIZED MASTER EQUATION 

We consider a diffusive process where particles are subject to annihilation using the microscopic approach of 
Einstein's original random walk model. For simplicity we consider a one-dimensional lattice where the particle hops 
to the nearest neighboring site (left or right) in one time step, and can then also be annihilated by a some reactive 
process as described by the discrete equation 

n*(r; t + 1) = ri*(r + 1; t _) + C+ 7i*(r - 1; t _) - n*(r; t +) , (3) 

where the Boolean variable n*(r; €) = {0, 1} denotes the occupation at time t of the site located at position r and ^± 
is a Boolean random variable controlling the particle jump between neighboring sites (^+ + ^ l)i while (^b, is the 
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reactive Boolean operator controling particle annihilation. The mean field description follows by ensemble averaging 
Eg.® with {n*{r;t)) = n{r;t), {£,±) = Pi, and {^n) = Ri, where i is an index for the position; using statistical 
independence of the ^'s and n*, and extending the possible jump steps over the whole lattice, we obtain 

+00 

n{r;t + dt) = Pj{r ~ j6r;t)n{r ~ jSr;t) - Rir;t) n{r;t) , (4) 

j = -oo 

where Pj (r — j6r) denotes the probability of a jump of j sites from site r — jSr, and R (r) the annihilation probability 
at site r; the number density is n (r; t) so that n (r; t) dr is the expected number of particles to find in the interval 
[r — dr/2,r + dr/2]. Note that in a closed system, i.e. without the second term on the right, the total number of 
particles, N, is constant so that one can divide through by this number to express the master equation in terms of 
/(r, <) = n(r,t)/N , the probability density. Alternatively, if the system contains multiple components, then a more 
useful concept is the concentration. For example, if there are two components, one of which is the solvent and the 
other the solute, then the solute concentration would be c(r, i) — n(r,t)/ (n{r,t) + ns{r,t)), where Ugir^t) is the local 
number density for the solute. In the common case that the solvent is uniform and stationary, ns{r,t) = Us, and that 
the solute is relatively dilute, rig >> n(r,t), one has, to first approximation, c{r,t) = n{r,t)/ns, which is what we will 
use in the following. 

In the classical case, the jump probabilities are constants, Pj{r — jdr;t) = Pj > with J^JL-ooPj = 1j is the 
reaction probability, R {r; t) = pa with 1 > Pr > 0- We take into account the configurational complexity of the 
reactive medium by allowing for the possibility that both the jump probabilities and the reaction probability are 
modified by interaction between the particles. This is modeled by writing Pj (r — jSr; t) = pj F [c{r — jSr; t)], with 
j 7^ and R (r; t) = prG [c (r; t)] giving the Generalized Master Equation 

+00 

c{r-t + 5t) -c{r-t)^ ^ {p^F [c{r ~ j5r;t)]{c{r ~ jSr-t) - c{r-,t))) - pRG[c{r-t)] c{r-t) . (5) 

j = -QO 

Notice that in order to retain their nature as probabilities, the functions F [c] and G [c] must both be greater than 
zero and less than one for all values of their arguments. 



III. DIFFUSION AND REACTION 

A. Generalized diffusion equation 

Considering the diffusive process alone, it was shown Q that the generalized diffusion equation that follows from 
Eq.(l5]) (without the second term on the r.h.s) is 



dt dr dr \ dx dy J ^ &■ 



Jr 



C^6t d ( dxF{x,y) dxF {x,y) fdxF{x,x)V\ dc 
2 dr \ dx dy \ dx J I dr ' 



with the compact notation (...)^ = {■■■)x=c{r t) y=c(r t)- Here C = [J^jjPjj IF advection speed and 

D = {j2j j^Pj ) ^^2St is diffusion coefficient. In Q it was also shown that the existence of a scaling solu- 
tion c (r; t) — t~^^^(l) {''"/t^^^) demands that F[c] ^ c^ in which case the scaling exponent is 7 = 2T1) ! since the jump 
probabilities Pj = pjF [c] must be < 1, one must have 77 > 0, that is 7 < 1, which is the signature of sub-diffusion 0. 
We now combine the description of sub-diffusion (with no drift, i.e. C = in ([6])) with reactive processes. 



B. Scaling reaction-diffusion 

Starting from the generalized master equation ([S]), we proceed along the lines of derivation of the generalized 
diffusion equation given in Performing a multiple scale expansion up to second order, we obtain the general form 



^ For j; = 0, one has 7 = 1, i.e. classical diffusion. The case of super-diffusion will be presented elsewhere. 



4 



of the reaction-diffusion (R-D) equation (with no drift and with reaction rate k = pji ^): 

-c (r; t) = (F[c (r; i)]c (r; t)) - k G[c (r; t)] c (r; t) . (7) 

As for the gcnerahzcd diffusion equation Q, we ask under which conditions there is a scaUng solution to equation 
([7]) of the form c (r; t) — t^''/^^ (r/f/^) = t^'^/^tj) (x). Expressing the time and space derivatives in terms of x, Eq.([7]) 
can be written as 



- l-^x<f> {x) = 2 Dt^'^^F (t-^''^(t> (x)) (j){x)- ktG (x)) (/. (x) 

The time-dependence on the right can only be eliminated if F{c) and G(c) have a functional power law form: F{c) 



(8) 



1 = ^1-^^(1-") 7/2^ and 1 = t ''/^ that is 



^(l-n)7/2 0n-l^ for some numbers a > 1 and n > 1 ; hence we must have 
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a + l 



n — 1 = — 

7 



(9) 



Thus, according to scaling consistency the exponents should be such that n = a + 2. When a > 1, we have anomalous 
diffusion: (r^) ^ t^+r (and more generally (r™) ~ and the reaction term goes like ~ —kcp". More exphcitly, 

using in ^ the reduced variable 



we obtain the scaled equation 



1 



:(C0(C)) - 0"(C)=O, 



(10) 



(11) 



(12) 



a + 1 dC 

which can be rewritten in terms of the original variables (r and t) to give 

Without the reactive term, i.e. with fc = 0, this reduces to our previous generalized diffusion equation in the absence 
of drift. Eauation (|12p is the generalized reaction-diffusion equation. 

IV. STEADY-STATE DISTRIBUTIONS 

In this Section we explore Eq.()12p as a natural extension of our previous description of generalized diffusion to 
include extinction. Because we are not solely interested in scaling solutions in this case, we will allow for arbitrary 
exponents a > 1 and n > 0. 



A. Boundary conditions 



One frequently studied problem is that of a semi-infinite sytem with constant injection of particles at the boundary. 
To be specific, we use the interval [0, oo] and note that the rate of change of the total number of particles in is simply 



dt 



Dn 



dt 
dc" (r; t) 



dr 



- Dn. 



dc" (r; t) 



dr 



- kus 



(9c" (r; t) 
di 



dr . 



r=0 "'O 

The first term on the right is the rate at which matter leaves the system via the boundary at infinity: we will assume 
that the concentration goes to zero sufficiently fast at infinity so that this term is zero - an assumption that will have 
to be checked a posteriori. The second term on the right is the rate at which particles are injected at the left boundary 
and the last term is the rate at which particles are removed by the extinction process. Our boundary condition will 
be to control the rate at which particles are injected so we set 



dN (t) 
dt 



= jo = -Dn, 



ac" (r; t) 



dr 



(14) 



r=0 



as the boundary condition of interest. 
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B. Steady-state solution 



We now seek a steady state solution with this boundary condition, 



= Z?|^c" (r) - fcc" (r) 



with 



9c" (r) 



dr 



r=0 



It is convenient to rewrite the problem with the change of variables 



Jo Jo = 



30 



so that the steady state equation has the simple form 



dz 



-^g{z) =g^ (z) with 



/kD 



dgiz) 



g ^ c 



dz 



-Jo- 



2 = 



This is integrated to get 



dgjz) 
dz 



±xlA+^g^iz). 
" a + n 



(15) 



(16) 



(17) 



(18) 



Recall that we assumed that the flux at infinity goes to zero. This means that either A = and limz^oo g (z) ~ or 
that A < and lim2_j.oo g (z) is finite . We rule out the latter case on the ground that without extinction we should 
get purely diffusive behavior and that adding extinction should not cause an increase in particles far from the source. 
A second integration then gives the implicit solution 



± z 



dg 



J^g' 



a + n 2a 
2a a — n 



5°2c" (z) -gV" (0) 



or, upon rearrangement, 



(19) 



.(.) = .(o) i±.--(o)^y^. 



(20) 



The boundary condition is 



.70 



dz 



2 a o+Ti , ^ 

—— g— , 

a + n 



Since we are interested in the circumstance that the injection rate is positive we must take the lower sign so 

' a + n\ 



g{z)=\ro 



a + n 
2a 



1 - 



zo 



, Zo = ajo 



or, rewriting the result in terms of the physical variables. 



c r 



Jo 



, a + n 



2a 



a~ n r 



2 ro 



, »'o = ajo 



2a 



a + ri,\"+" ID 
2a ) VT 



(21) 



(22) 



(23) 



There are two cases that must be distinguished depending on whether n > a or a > n. In the first case the solution 
has infinite support and is a simple algebraic decay 



c{r)= Jo 



a + n 
2a 



n — a r 



2 ro 



n > a . 



(24) 
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The second, more complicated case occurs when a > n. Then it is clear from Eg. ([23]) that the concentration will, in 
general become imaginary and in all cases its magnitude will increase without bound for sufficiently large r. The only 
way to avoid this unphysical behavior is if the solution has finite support so that 

where the step function (x) = 1 for a; > and zero otherwise. Noting that 

^fix)e{x)^f'{x)Q{x) + f{0)d{x) , 
ax 

(x) e (x) ^ /" (x) e (x) + r (0) Six)+f (0) 6' (x) , 

it is clear that (|25p can only be an acceptable solution to the steady state equation ([T5|) if the first two derivatives 
of the coefficient of the step function vanish at r = . This simply imposes the requirement on the exponent that 

— 2 > which is always true provided that ti > 0, as was already required. Thus, the final, physically valid 
solution with finite support ([55]) is restricted to a range of values of the coefficients a > n > 0. 
We note that the solutions (|24l) and ([25)) can be expressed as g-exponentials, 

Cq (x) — (1 + (1 — g)a;)^"'' 8(1 + (1 — f?)^;) with the identification q = _|_ i and that q > 1 gives the case of 
infinite support while 1 > q gives the case of finite support. From the properties of the (/-exponential we know that for 

g = 1 the decay of the concentration will be expontial, c (r) = (jg) ° e~^'^° with tq = a^i This of course includes 
the steady state solution of the classical reaction-diffusion equation with a = n = 1. 

The physical interpretation of these results can be understood as follows: increasing n decreases the extinction rate 
(since the reaction term goes like c" and c < 1) while increasing a decreases the rate of diffusion (this is easily seen 
from the scaling r ^ t'^^'^ or by writing the diffusion term as -^D^^ = ^ (aDc"^^) so the effective diffusion 
coefficient goes like c"~^). Hence, making n large or a small leads to infinite support: diffusion is fast, extinction is 
slow. The converse, making n small or a large leads to finite support because diffusion is slow and extinction is fast. 
The resulting steady state profiles are compared in Fig. [TJ 



C. Robustness of the steady state 

The question of robustness is an important issue as discussed by Eldar et al. [§| and by Yuste et al. Q in particular 
for morphogen gradient formation as precursor to cell differentiation. Robustness is a measure of the strength of the 
steady state profile versus changes in the variables controlling input flux and degradation, such as jo and k. The cited 
authors characterized it as the quantity TZ}, = d \dL/d\ogb\~^ where d is a characteristic microscopic length (e.g. the 
cell size) and b denotes jo or k; L is the distance at which the steady state c(r) takes a given value and is obtained by 
inversion of the steady state solution c{r)r=L- A high value of TZb is an indication of the buffering capacity against 
changes in the input flux and degradation rate. Here, however, we prefer to consider directly the relative change in 
the concentration at point r due to a change in the value of quantity 6, thereby dcflning the (position-dependent) 
sensitivity to parameter b as 

SUr) = '-^. (26) 
For n > a, the case of infinite support, a short calculation gives the sensitivity as 

'^.oW-^^YT^ ; n>a, (27) 
and for n = a = 1, i.c in the classical case of exponential decay, this becomes 

5jo(r) = l ; n = a = l, (28) 

which we will take as a reference point. One also gets exponential decay for the more general condition n = a (see 
section BVBp . but in this case we find 

5jjr) = i ; n = a, (29) 
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so that even though the decay is exponential, it is nevertheless true that increasing the nonlinearity of the process 
decreases the sensitivity of the concentration to variations in the injection rate. Note that the general result for 
infinite support is bounded by 

5,0 (r)<^— ; n>a, (30) 
a + n 

so that - independent of position - increasing nonlinearity in either the diffusion process or in the extinction process 
has the effect of buffering the concentration against changes in the rate at which material is injected. 
The case of finite support, a < n, is more complicated. A simple calculation gives 

5.0 M = Y-^^^ ; a>n, (31) 

2 rg 

SO that there are two effects at work: decreasing sensitivity with increasing nonlinearity, as above, and increasing 
sensitivity with increasing distance from the source. In fact, in this case we find 

(r) > 1 ^ r > = (l —) tq. (32) 

a — n \ a + n J 

Clearly, this is only relevant if the right hand side is less than tq. For n < 1, this is always the case: i.e., there is 
always a region of enhanced sensitivity in the range < r < vq. For n > 1, there is a region of enhanced sensitivity 
for 



a > a* = 1 + J(n - 1) (n + 3) = n + 2 h ... ; a > n > 1. (33) 

n 

Only for the restricted range a* > a > n > 1 is there no region of enhanced sensitivity for the case of finite support. 

In summary, we find that (i) for infinite support, n > a, increasing nonlinearity always decreases sensitivity of the 
concentration to the injection rate; (ii) the same holds true for the case of finite support when a^, > a > n > I; (iii) 
the case of finite support will, for n < 1 or a > a* show enhanced sensitivity in the region < r < ro . 



V. COMPARISON TO SIMULATION AND EXPERIMENT 



A. Numerical solution of master equation 

We have performed numerical computation of the master equation ([5]) in order to verify three aspects of this theory: 
first, that the non-linear dynamics eventually leads to a steady state; second, that the steady state is independent 
of the initial conditions and third, that our analytic, continuum result is a good representation of the steady state. 
Figure [5] shows the result of solving the master equation with an initial condition c(r) = and with constant flux 
at the origin for two cases: one with finite support, n < a, and one with infinite support, n > a. In both cases, we 
do indeed find that at long times the system settles into a steady state that is well-described by the analytic results, 
p4|) and ([25]). Note that, in the case of infinite support, one must go to somewhat longer times to reach the steady 
state. To test that the sensitivity of the steady state to the boundary conditions, the calculations were repeated with 
a boundary condition of fixed value of the concentration at r = 0. The result for the case of finite support is shown in 
Fig. |3]whcre it is again seen that the system reaches a steady state and that the steady state is that of the continuum 
theory. Similar results were found for the case of infinite support. This comparison of numerical and analytical results 
therefore shows good agreement between the continuum approximation and the discrete microscopic dynamics and 
furthermore provides evidence that the steady state is unique. 



B. Comparison to experiment 

As an application of the theory we compare our analytical solution for the steady state with experimental results 
obtained from measurements performed in the Drosophila wing disc where morphogens are produced by a subset of 
cells wherefrom they diffuse and are degraded thereby forming a concentration gradient whose profile shape appears 
crucial for subsequent cell specification This situation corresponds to the reaction-diffusion theory presented in 
the present article. Experimental results given in [loj present the intensity signal of the Wg morphogen as a function 
of distance from the source obtained by image processing showing the profile of the diffusing protein in selected regions 
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FIG. 1. Steady state: c(r)/c(0) = exp(-r/ro) (g = 1; black); c(r)/c(0) = + ^^7^,) " ° for n = a + 2 (infinite support, 
g = 2; red) and a = ?i + 0.8 (finite support, q — 0.6; blue). 



0.2. 




r/8r 



FIG. 2. Numerical and analytical solutions for the steady state profile. Left panel: Case of finite support. Numerical solution 
of the Master Equation ((5| where pj F = pj c"~^ with pj = 0.2 for j G [—2, +2] and a = 1.5 and G = pR c" with n — 1 and 
PR = 10"^, for t = 50, 200, 500, 800, 1.5 x 10^, 2 x 10^, 1.5 x lO" time steps (symbols); the boundary condition is finite fiux at 
r = and the initial condition is zero concentration everywhere. Analytical steady state solution (|25|l (black curve). Right panel 
: Case of infinite support. Same as left panel except n = 2 and pR = 10^^ for t = 50, 200, 500, 10^, 2 x 10^, 4 x 10^, 10"*, 1 x lO'^ 
time steps (symbols) and steady state solution (|24|) (black curve). Note that there are no adjustable parameters in either case. 

of the Drosophila wing disc. In the absence of numerical data, we processed the signal images to obtain the data 
shown in Figs. S] and [5] where they are compared to our analytical results. Clearly we find that the sub-diffusive 
nonlinear reactive steady state profile (pi)) with infinite support reproduces very well the experimental data indicating 
slow degradation combined with extended sub-diffusion. In all cases, we also show best-fits to an exponential of the 
form f{z) = Ae^~^l^l^ and it is clear that the experimental data are very poorly fit by an exponential decay. 

VI. COMMENTS 

Wc derived the nonlinear reaction-diffusion equation starting from Einstein's microscopic model where the diffusing 
particles are also subject to an annihilation reactive process. The nonlinear reaction-diffusion equation was obtained 
under the demand that scaling be satisfied for diffusive motion wherefrom a relation follows between the scaling 
exponent and the nonlinear exponents whose range of possible values exhibit the signature of sub-diffusion. While full 
scaling should in principle be satisfied for the space-time dependent equation, this requirement can be relaxed between 
the reaction term exponent and the scaling exponent for the steady state equation. This observation is important for 
the analysis of the R-D steady state solutions which take the form of a power law with in one case infinite support 
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FIG. 3. Numerical and analytical solutions for finite support profile (n < a) with boundary condition of fixed c(0). Left panel: 
Numerical solution of the Master Equation ([5]) where pj F = pj c"~^ with pj = 0.2 for j = [—2, +2] and a = 1.5 and G — pr c" 
with n = 1 and pfl = 10"^, for t = 3 x 10^,6 x 10^,8 x 10^,1.2 x 10^,6 x 10^ time steps. Right panel: Comparison between 
numerical solution of the Master Equation (O for t = Qx 10^ time steps (open circles) and analytical steady state solution (|25p 
(black curve). Note that there are no adjustable parameters. 




1.5 -0.25 0.25 0.5 

Distance 



FIG. 4. Experimental data (black dots) from Han et al, Fig. 6. A in [T^ from the fluorescence intensity of the Wg protein 
(vertical axis; normalized values) versus distance (in a.u.; horizontal axis) measured from the anterior-posterior axis along the 
dorsoventral direction in the posterior compartment of the Drosophila wild- type wing disc • The black curve is the best-fit 
of the theoretical steady state with n — a ~ 3.8. For comparison the dashed curve shows the best-fit exponential profile. 



and in the other case finite support. 

We discussed the sensitivity of the steady state versus changes in the input flux and we found that profiles with 
infinite support show minimal sensitivity, and such profiles with infinite support were shown to correspond to experi- 
mental observations. On the other hand we showed that profiles with finite support should exhibit stronger sensitivity 
to input fiux changes, and it seems that such profiles with finite support have not been observed in morphogen gradient 
formation. This observation may suggest that extreme sensitivity excludes this type of profile in natural morphogen 
gradient formation because degradation is too fast with respect to diffusion in order to establish the necessary gradient 
for subsequent cell differentiation. 
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FIG. 5. Same as Fig|3]for experimental data (black dots) from Han et al, Fig.6.B in [T3| for a mutant strain. Fit of theoretical 
steady state (|24p (black curve) to the experimental data; because of the obvious asymmetry of the data along the dorsoventral 
axis, the left panel shows a fit based only on the data for negative distances, giving 71 — a ~ 1.7, and the right panel shows a 
fit to data for positive distances, giving n — a ~ 3.3. In both cases, a best-fit to an exponential decay is shown as the dashed 
curves. 
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